#load data and fit models
 
tadpoles <- read.table("ecol 563/tadpoles.csv", header=T, sep=',')
tadpoles$fac3 <- factor(tadpoles$fac3)
out1 <- lm(response~fac1 + fac2 + fac3 + fac1:fac2 + fac1:fac3 + fac2:fac3 + fac1:fac2:fac3, data=tadpoles)
out2 <- lm(response~fac1 + fac2 + fac3 + fac1:fac2 + fac2:fac3, data=tadpoles)
myvec <- function(fac1, fac2, fac3) c(1, fac1=='No', fac1=='Ru', fac2=='S', fac3==2, (fac1=='No')*(fac2=='S'), (fac1=='Ru')*(fac2=='S'), (fac2=='S')*(fac3==2))
 
# calculate treatment means using model
fac.vals2 <- expand.grid(fac1=c('Co','No','Ru'), fac2=c('D','S'), fac3=1:2)
cmat <- t(apply(fac.vals2, 1, function(x) myvec(x[1],x[2],x[3])))
ests <- cmat%*%coef(out2)
vmat <- cmat %*% vcov(out2) %*% t(cmat)
se <- sqrt(diag(vmat))
low95 <- ests + qt(.025,out2$df.residual)*se
up95 <- ests + qt(.975,out2$df.residual)*se
results.mine <- data.frame(fac.vals2, ests, se, low95, up95)
results.mine$fac1.num <- as.numeric(results.mine$fac1)
 
# generate graph from last time for sibship 1
plot(c(.8,3.2), range(results.mine[,c("up95", "low95")]), type='n', xlab="Hormonal treatment", ylab='Mean mitotic activity', axes=F)
axis(1, at=1:3, labels=levels(results.mine$fac1))
axis(2)
box()
results.mine.sub1 <- results.mine[results.mine$fac3==1,]
#first food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='D',]
#amount of shift
myjitter <- -.05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=1)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=1)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=1, pch=16)
#second food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='S',]
myjitter <- .05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=2)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col='white', pch=16, cex=1.1)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2, pch=1)
mtext(side=3, line=.5, 'Sibship 1', cex=.9)
legend('topleft',c('Detritus','Shrimp'), col=1:2, pch=c(16,1), cex=.8, pt.cex=.9, title=expression(bold('Food Type')), bty='n')
summary(out2)
 
 
# function to calculate difference-adjusted confidence intervals
ci.func <- function(rowvals, lm.model, glm.vmat) {
nor.func1a <- function(alpha, model, sig) 1-pt(-qt(1-alpha/2, model$df.residual) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*%c (1,-1)), model$df.residual) - pt(qt(1-alpha/2, model$df.residual) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*% c(1,-1)), model$df.residual, lower.tail=F)
nor.func2 <- function(a,model,sigma) nor.func1a(a,model,sigma)-.95
n <- length(rowvals)
xvec1b <- numeric(n*(n-1)/2)
vmat <- glm.vmat[rowvals,rowvals]
ind <- 1
for(i in 1:(n-1)) {
for(j in (i+1):n){
sig <- vmat[c(i,j), c(i,j)]
#solve for alpha
xvec1b[ind] <- uniroot(function(x) nor.func2(x, lm.model, sig), c(.001,.999))$root
ind <- ind+1
}}
1-xvec1b
}
 
# confidence levels for all 12*11/2 = 66 mean comparisons
ci.func(1:12, out2, vmat)
 
# look only at first six means, the ones for sibship 1
fac.vals2
ci.func(1:6, out2, vmat[1:6,1:6])
 
# calculate 75% and 83.5% confidence intervals
low95 <- ests + qt(.025,out2$df.residual)*se
low75 <- ests + qt((1-.75)/2,out2$df.residual)*se
up75 <- ests + qt(1-(1-.75)/2,out2$df.residual)*se
low83 <- ests + qt((1-.835)/2,out2$df.residual)*se
up83 <- ests + qt(1-(1-.835)/2,out2$df.residual)*se
 
# add them to the data frame
results.mine <- cbind(results.mine,low75,up75,low83,up83)
results.mine
 
# add difference-adjusted cis to the graph
par(lend=1)
plot(c(.8,3.2), range(results.mine[,c("up95", "low95")]), type='n', xlab="Hormonal treatment", ylab='Mean mitotic activity', axes=F)
axis(1, at=1:3, labels=levels(results.mine$fac1))
axis(2)
box()
results.mine.sub1 <- results.mine[results.mine$fac3==1,]
#first food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='D',]
#amount of shift
myjitter <- -.05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=2)

#difference-adjusted confidence intervals
segments(cur.dat$fac1.num+myjitter, cur.dat$low83, cur.dat$fac1.num+myjitter, cur.dat$up83, col='lightpink3', lwd=5)
segments(cur.dat$fac1.num+myjitter, cur.dat$low75, cur.dat$fac1.num+myjitter, cur.dat$up75, col='red4', lwd=8)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col='white', pch=15, cex=1.2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2, pch=22, cex=1.2)
#second food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='S',]
myjitter <- .05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=4)

#difference-adjusted confidence intervals
segments(cur.dat$fac1.num+myjitter, cur.dat$low83, cur.dat$fac1.num+myjitter, cur.dat$up83, col='lightblue', lwd=5)
segments(cur.dat$fac1.num+myjitter, cur.dat$low75, cur.dat$fac1.num+myjitter, cur.dat$up75, col='dodgerblue4', lwd=8)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=4)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col='white', pch=16, cex=1.2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=4, pch=1, cex=1.2)
mtext(side=3, line=.5, 'Sibship 1', cex=.9)
legend('topleft',c('Detritus','Shrimp'), col=c(2,4), pch=c(22,1), cex=.8, pt.cex=.9, title=expression(bold('Food Type')), bty='n')
 
# the seq function generates a vector of number with equal increments
seq(3.3,4.5,.05)
# add grid lines to graph
abline(h=seq(3.3,4.5,.05), col='grey70', lty=3)
 
# repeat graph with only one set of difference adjusted cis
par(lend=1)
plot(c(.8,3.2), range(results.mine[,c("up95", "low95")]), type='n', xlab="Hormonal treatment", ylab='Mean mitotic activity', axes=F)
axis(1, at=1:3, labels=levels(results.mine$fac1))
axis(2)
box()
results.mine.sub1 <- results.mine[results.mine$fac3==1,]
#first food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='D',]
#amount of shift
myjitter <- -.05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=2)

#difference-adjusted confidence intervals
segments(cur.dat$fac1.num+myjitter, cur.dat$low83, cur.dat$fac1.num+myjitter, cur.dat$up83, col='lightpink3', lwd=5)
#segments(cur.dat$fac1.num+myjitter, cur.dat$low75, cur.dat$fac1.num+myjitter, cur.dat$up75, col='red4', lwd=8)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col='white', pch=15, cex=1.2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=2, pch=22, cex=1.2)
#second food group
cur.dat <- results.mine.sub1[results.mine.sub1$fac2=='S',]
myjitter <- .05
arrows(cur.dat$fac1.num+myjitter, cur.dat$low95, cur.dat$fac1.num+myjitter, cur.dat$up95, angle=90, code=3, length=.05, col=4)

#difference-adjusted confidence intervals
segments(cur.dat$fac1.num+myjitter, cur.dat$low83, cur.dat$fac1.num+myjitter, cur.dat$up83, col='lightblue', lwd=5)
#segments(cur.dat$fac1.num+myjitter, cur.dat$low75, cur.dat$fac1.num+myjitter, cur.dat$up75, col='dodgerblue4', lwd=8)
lines(cur.dat$fac1.num+myjitter, cur.dat$ests, col=4)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col='white', pch=16, cex=1.2)
points(cur.dat$fac1.num+myjitter, cur.dat$ests, col=4, pch=1, cex=1.2)
mtext(side=3, line=.5, 'Sibship 1', cex=.9)
legend('topleft',c('Detritus','Shrimp'), col=c(2,4), pch=c(22,1), cex=.8, pt.cex=.9, title=expression(bold('Food Type')), bty='n')
abline(h=seq(3.3,4.5,.05), col='grey70', lty=3)
 
# confidence levels for the second set of means
ci.func(1:6,out2,vmat[7:12,7:12])

#all pairwise comparisons
 
# function taken from web to create a matrix of 0s and 1s
p.table = function(x, g, p.adjust.method = "none", ..., level = 0.05) { 
  ## fill out p-value table 
  p = pairwise.t.test(x, g, p.adjust.method, ...)$p.value 
  p[is.na(p)] = 0 
  p = rbind(0, cbind(p, 0)) 
  p = p + t(p) 
  ## 0 = no difference, 1 = difference 
  p = 1 * (p <= level) 
  diag(p) = 0 
  ## get means and find their order 
  m = tapply(x, g, mean) 
  o = order(m) 
  p = p[o, o] 
  dimnames(p) = list(names(m[o]), names(m[o])) 
  cbind(mean = m[o], p) 
  }

#use function on treatment means
p.table(tadpoles$response,tadpoles$treatment)
 
# we need to remove missing values first
tadpoles2<-tadpoles[!is.na(tadpoles$response),]
dim(tadpoles)
dim(tadpoles2)
p.table(tadpoles2$response, tadpoles2$treatment)
 
#classify means into groups based on pairwise diffs table
diffs <- c('a','a','a','b','c','c','cd','cd','cd','de','ef','f')
 
#obtaining treatment means
out0 <- lm(response~treatment, data=tadpoles)
summary(out0)
out0a <- lm(response~treatment-1, data=tadpoles)
summary(out0a)
out0b <- lm(response~fac1:fac2:fac3-1, data=tadpoles)
summary(out0b)
 
#dynamite graph

# sort means in ascending order
my.means <- sort(coef(out0a))
my.means
 
# sort std errs in same order as the means
se.temp <- summary(out0a)$coefficients[,2]
se.temp2 <- se.temp[order(coef(out0a))]

# rescale means so error bars show up on graph
dat1 <- my.means-3
# create bar plot and save coordinates of bars
bp <- barplot(dat1, axes=F, ylim=c(0,1.6), names.arg=substr(names(my.means),10,nchar(names(my.means))),cex.names=0.8)
# add correct labels on y-axis (because of rescaling)
axis(2, at=seq(0,1.4,.2), labels=seq(0,1.4,.2)+3)
# add handles to dynamite plot
arrows(x0=bp,y0=dat1,x1=bp,y1=dat1+se.temp2,angle=90,length=0.1)
# add labels to the handles
text(bp, dat1+se.temp2, diffs, pos=3, cex=.9)